Astronomy & Astrophysics manuscript no. paper 
March 19, 2012 



©ESO2012 



Compression Behaviour of Porous Dust Agglomerates 



A. Seizinger/ R. Speith,^ and W. Kley' 



(N 

o 

(N 



^ 



6 



Institut fiir Astronomie and Astrophysik, Eberhard Karls Universitat Tubingen, 
Auf der Morgenstelle 10c, D-72076 Tubingen, Germany 
e-mail: alexsOtat . physik . uni -tuebingen . de 
Physikalisches Institut, Eberhard Karls Universitat Tubingen, 
Auf der Morgenstelle 14, D-72076 Tubingen, Germany 



Received 20.01.2012; accepted 11.03.2012 



ABSTRACT 



Context. The early planetesimal growth proceeds through a sequence of sticking collisions of dust agglomerates. Very uncertain is 

still the relative velocity regime in which growth rather than destruction can take place. The outcome of a collision depends on the 

bulk properties of the porous dust agglomerates. 

Aims. Continuum models of dust agglomerates require a set of material parameters that are often difficult to obtain from laboratory 

experiments. Here, we aim at determining those parameters from ab-initio molecular dynamics simulations. Our goal is to improve 

on the existing model that describe the interaction of individual monomers. 

Methods. We use a molecular dynamics approach featuring a detailed micro-physical model of the interaction of spherical grains. 

The model includes normal forces, rolling, twisting and sliding between the dust grains. We present a new treatment of wall-particle 

interaction that allows us to perform customized simulations that directly correspond to laboratory experiments. 

Results. We find that the existing interaction model by Dominik & Tielens leads to a too soft compressive strength behavior for uni- 

and omni-directional compression. Upon making the rolling and sliding coefficients stiffer we find excellent agreement in both cases. 

Additionally, we find that the compressive strength curve depends on the velocity with which the sample is compressed. 

Conclusions. The modified interaction strengths between two individual dust grains will lead to a different behaviour of the whole 

dust agglomerate. This will influences the sticking probabilities and hence the growth of planetesimals. The new parameter set might 

possibly lead to an enhanced sticking as more energy can be stored in the system before breakup. 

Key words. Planets and satellites: formation ~ Protoplanetary disks - Methods: numerical 
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1. Introduction 



Unraveling the question of planetesimal formation is a crucial 
issue for the core acc r etion model of planet formation pro- 
posed by iPollack et alJ (1199 6). To this day it remains unclear 
how dust and ice particles can grow several orders in magni- 
tude from micron to kilometer sized objects. In the beginning 
micron sized dust grains (further referred to also as monomers 
or particles) may grow by low velocity, hit-and-stick colli- 
sions resulting in h ighly porous fractal aggregates ( Kem pf et aU 
Il999t iBIum et all 200 0). As the aggregates grow larger their 
motion increasingly decouples fr om the surround i ng ga s lead- 
ing to higher collision velocities dWe idenschill ind, I l977h . With 
increasing impact energy colhding a ggregate s get cornpacted 
(Blum & Wurm, 2000; Suvamaetal.1 l2008t IWada et all l200i 
[Paszun & Dominik , .2009.) . At some point relative velocities be- 
come large_enoughtha^Jragmentation is supposed to set in 
(e.g. iBIum & WurmL |2008), which may limit the collisional 
growth of aggr egates. Apart from fragmentation, radial dr ift and 
bouncing (e. g. Langkowski et all l2008t IWeidling et al.L l2009t 
iGuttler et al.L 20101) may hamper the growth of planetesimals 
dZsom et al.L l2010t) . Despite these obstacles, experimental evi- 
dence indicates possible growth in high speed collisions (a few 
lOm /s) through sticking and reaccretion of ejecta ( Wurr net al.L 



12005 ; Teiser & Wurm, 2009) or the sweep-up of smaller parti- 
cles (iWindmarket al.i.i2012) . To understand the possible growth 



regimes requires material properties and simulations beyond the 
current data base. 

In the context of planetesimal formation, collisions of mi- 
cron sized aggregates have been st udied theoretically using a 
molecular dynam ics approach (e.g. iDominik & Tieletii . 119971 ; 
IWada et alll2007l) . Here, the motions of all grains that make up 
the aggregate are followed individually, considering suitable in- 
teraction forces. However, billions of such grains would be nec- 
essary to model meter-sized objects. For computational reasons 
it is therefore necessary in such cases to use continuum models 
to simulate collisions between macroscopic aggregates. Smooth 
Particle Hydrodynamics (SPH) constitutes such an approac h 
dSirond l2004t ISchafer et all 120071; iGeretshauser et all l2010l) . 
However, the outcome of such simulation s depends s trongly on 
the underlying porosity model (Guttler e t al.L l2009h . To sim- 
ulate collisions of large dust boulders with SPH, several pa- 
rameters describing the behaviour of the material such as the 
compressive-, tensile-, and shear-strength must be known in ad- 
vance. 

To this day, only few attempts have been made to measure the 
required material parameter of porou s dust aggregates in labora- 
tory experiments. , Blum & Schrapleij d2004 ) studied the case of 
unidirectional compression, where a sample aggregate is com- 
pacted between two plates. While the bottom plate was fixed a 
certain load was applied to the top plate. To calculate the pres- 
sure the cross-section of the compressed aggregate had to be de- 
termined. Their samples were produced by random ballistic de- 
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position (RBD) and featured an initial filling factor of 0.15. The 
maximum filling factor they could reach was limited to approx- 
imately 0.33 as grains began to flow outwards as the pressure in 
the center inc reased. 

Later on, iGiittler et alj (l2009l) used a similar approach to 
study the omnidirectional compression. Their sample was put 
into a solid box and the load was applied by a movable piston 
on the top (see Fig.|3]l. As the grains could not escape the box 
they reached a significantly higher filling factor of roughly 0.58 
upon compression. One advantage of this approach lies in the 
elimination of any uncertainty in the determination of the fill- 
ing factor since the volume currently occupied by the sample is 
unambiguous. 

The compressive stren gth of porous dust aggreg ates was 
determined numerically by iPaszun & DominiM ( l2008l) using a 
mo lecular dynamics approach based on an interaction model 
by iDominik & TielensI (I1995L [1996). Yet, they only modeled 
the unidirectional compression and their simulations were lim- 
ited to very few particles {^ 300). Using similar material 
and model parameters the model has been applied to explore 
the growth regim es of dust a gglornerates unde r protoplane- 
tary conditio ns (Suva maetall EoOSL IWada et all l2009l boiU 



Table 1. Material Parameters. 



IOrmeletal.L[2009) . 

In this work we step back again and present a method to 
obtain the continuum material parameter of porous particle ag- 
glomerates, in particular the compressive strength curve, from 
ab-initio simulations using a molecul ar dynamics approach. Our 
appr oach is based on the model bv IDominik & TielensI (1 19951 
11996 ) using extensions by Wad a et aLl (l2007h . To test the ap- 
plicability of the model we perform customized simulations of 
both, omni- and unidirectional compression, with a much greater 
number of particles of the order of lO'*. As we will see, modifi- 
cations of the model ar e required to properl y reproduce the ex- 
perimental results from IGiittler et al.| (l2009l) . Having calibrated 
our model to the case of the slow compression as measured by 
IGiittler et al.l 02009) we will subsequently present new results on 
how the compressive behaviour changes with increasing com- 
pression speed. 

First, we briefly summarize the underlying physical model 
and present our extensions in Sect. 2. Our numerical setup is ex- 
plained in Sect. 3. In Sect. 5, we first describe the calibration of 
our model. Afterwards, we present our results of studying the 
dynamic compression and provide simple analytical approxima- 
tions describing the dependence of the compressive strength on 
the compression speed. 



2. Interaction model 

The agglomerates used in our simulations are composed of 
spherical monomers of equal size. To describe the interaction 
between individual monomers we adopt the physical model 
that has been presented by Dominik & Tielens (1997). In this 
model elastic grains may establish adhesive contacts when 
touching each other. Upon deformation of these contacts ki- 
netic energy is dissipated. The forces and torques acting upon 
the mono mers can be derived from coiTesponding potentials 
dWada et al., 2007i) . 

2.1. Particle-Particle interaction 



In lDominik & TielensI d 19971) the interaction between two spher- 
ical particles is divided into four types (see Fig.lTjl: 

1 . Compression/ Adhesion 



Pliysical property 



Silicate 



Particle Radius r (in yum) 

Density p (in g cm"^) 

Surface Energy y (in mj m"^) 

Young's Modulus E (in GPa) 

Poisson Number v 

Critical Rolling Length fcrit (in nm) 



0.6 

2.65 

20 

54 

0.17 

2 



2. Rofling 

3. Sliding 

4. Twisting 



In accordance with IWada et aP (l2007l) the following notation 
will be used in this work: r, denotes the radius of a monomer 
/, y, the surface energy per area, £, the Young's modulus, y, the 
Poisson's ratio, and G, the shear modulus. Furthermore, the re- 
duced radius R is given by 



and 



n + r: 



LP* — \ ' . ■' 



.2^1 



Gh 
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2 - V? 2-v 
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G, 



G, ) 



Note that for simplification all monomers in our simulations fea- 
ture the same properties. 

2.1.1. Material Parameters 

The basis of our simulations are the material parameters of sili- 
cate as summarized in Tab.[T] With the exception of the surfac e 
energy y these values comply with iPaszun & DominiH (l2008l) . 
where they used y = 25 mJ/m^. These parameters are also 
in reasonable agreement with experimental d ata as quoted by 
iBlum & Schrapleil (l2004 : IGiittler et al.1 (l2009l) . In this work we 
use a slightly lower val ue of y = 20 mJ/m^ w hich agrees with 
recent measurements of iGundlach et al.l (1201 lb . 

2.1.2. Compression/Adhesion 

The par ticle interaction in the normal direction has been devel- 
oped bv Johnson et al.l (Il97ll) (often refeiTed to as JKR theory). 
They extend the Hertzian theory by taking adhesion due to sur- 
face forces into account. 

To very briefly summarize their model let us first note that 
the compression length 6 of two monomers, which are in contact 
with each other and located at jci and x^, is defined by 



d^n +r2-\\x\ -Xt\\ 



(1) 



where ||u|| denotes the norm of the vector u. The radius a of the 
corresponding contact area can be obtained by 



5o \ao/ \ao 



1/2 



(2) 
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Fig. 1. The four types of particle interaction: Compression/ Adhesion (a), Rolling (b), Sliding (c), and Twisting (d). The left panel 
depicts particle-particle interactions while on the right the corresponding particle-wall version is displayed. 



where 6q and gq denote the equilibrium compression length and 
contact radius, respectively, where the repulsive and attractive 
normal forces are equal. It applies 



. "0 


' 






The force acting upon the monomers is 


F = 4Fc 


\flO/ \flo/ 


' 



(3) 



where Fc - inyR is the force needed to break the contact. 

A new contact is established if two freely moving particles 
touch each other, which means 6 - Q. However, once a contact 
has formed it will not break until the compression length ex- 
ceeds a certain threshold 6c - (9/16)'^"^ Sq. Thus, contacts may 
be stretched out a little bit before they finally break. 



2.1.3. Sticking velocity 

On contact formation or breaking there is a jump in the potential 
corresponding to the JKR force (see Wada et al., 2007), which 
reflects the dissipation of k i netic energy upon contact formation 
or breaking. IChokshi et al.l (Il993h proposed that this energy dis- 
sipation is caused by the excitation of elastic waves. From this 
energy dissipation one can calculate the maximum velocity Vstick 
at which particles stick on head-on impacts in JKR-theory 



Vstick = 1-07 



r 



5/6 



E*ViR5/6p 



1/2- 



(4) 



For micron-sized Si02 grains, we obtain velocitie s of the or- 
der o f 0.1 m/s. However, laboratory experiments by Poppe et alJ 
( 12000.) yield a significantly higher sticking velocity of roughly 
1.2ms ~' for similar sized partic les. To overcome this discrep- 
ancy, iPaszun & DominiM (l2008h introduced another damping 
mechanism that dissipated additional kinetic energy upon the 
contact creation and thus increased the sticking velocity to match 
the laboratory experiments. 

Since the sticking velocity has to be increased by approxi- 
mately one order of magnitude a significant amount of kinetic 
energy has to be dissipated aside from JKR-theory. Dissipating 
the kinetic impact energy all at once leads to a significant change 
of the relative velocity between the colliding particles. As they 



may be in contact with other particles as well, significantly mod- 
ifying the velocity of one particle during one integration step 
may introduce numerical hazards. F or low collision ve l ocities 
(the compression velocity used by Paszun & Dominila, l2008l 
was 0.05 ms"^) the additional damping is low and therefore this 
problem does not arise. However, since in this work we also want 
to study the dynamic compression behaviour at high velocities 
we choose a different approach to adjust the sticking velocity as 
explained in the next section. 

2.1.4. Normal oscillations 

The normal force may be both attractive or repulsive depend- 
ing on the current compression length 6. When two particles 
come too close to each other they are repelled, whereas they 
get attracted to each other due to adhesive surface forces while 
the contact is stretched out. This leads to oscillations of two 
monomers in the normal direction of a contact (from now on 
referred to as normal oscillations). In reality, one would expect 
tha t these oscilla t ions a re eventually damped away. For instance 
Bril liantov et alJ (l2007h proposed a viscoelastic damping mech- 
anism. 

From a numerical point of view these normal oscillations 
constitute a major nuisance. Not only may they artificially heat 
up aggregates (Paszun & Dorninik, 2008) but most importantly 
they need to be properly resolved in time. For micron-sized Si02 
grains the typical timescale of these oscillations is of the or- 
der of 10ns. Therefore our integration timestep is limited to 
=a 1 -0.3ns. 

iPaszun & DominikI (|2008|) tackle the problem of artificial 
heating by introducing an additional, weak damping mechanism 
which has hardly any influence on the sticking velocity but even- 
tually damps away normal oscillations. In this work we follow a 
similar approach. As before we consider two monomers located 
at Xi and X2 that are in contact with each other. The contact nor- 
mal «c is 

"': = t; u ■ (5) 

||jCi - X2\\ 

The relative velocity Vrei in normal direction of the contact is 
then given by 



Vid = (Vi - V2) • «c 



(6) 



where Vi, V2 denote the velocity of particle 1 and 2, respectively. 
The viscous damping force f damp is 

Fdamp = -A-v'ieinc , (7) 

where k is a damping constant determining the strength of the 
damping. We use /c = 1 x 10"* g/s. 
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2.1 .5. Rolling, sliding, and twisting 

The forces and torques resulting from the t angential motion of 
the contact area have been formulated by Dominik & Tielens 
( 1995, 1996). So called "contact pointers" (Dominik & Nubold, 
r2002h provide a convenient way to track the evolution of the 
contact area. They are unit vectors that initially point from the 
center of the particle to the center of the contact area. Due to 
the rotation of the particles their orientation changes over time 
(JDominik & Nubold, 2002). Contact pointers are used to define 
the rolling, sliding, and twisting displacement which quantify 
how much rolling, sliding, or twisting of a contact has occurred. 

Using these displacements, Wada et al. (2007) derived the 
forces and torques which agree with Dominik & TielensI (1 19951 
[1996), from corresponding potentials. All three types of inter- 
action remain elastically as long as the displacements do not 
exceed certain thresholds (from now on referred to as critical 
displacements). 

Let «i and «2 be the corresponding contact pointers describ- 
ing the contact. The rolling displacement^, sliding displacement 
^, and twisting displacement are then defined by: 



^ = /? (Ml + M2) , 

f = riiti - r2n2 - (fiiii ■ n^ - rifii ■ 11^^)11^ 
4> = «c(0 I (cJiiO - co2(t')) ■ n^(t')dt' , 

J to 



(8) 
(9) 

(10) 



where fo denotes the time when the contact has been established. 
For the rolling motion, the forces f , and torques Af, acting 
upon particle 1 due to being in contact with particle 2 are 

fr = 0, (11) 

Mr = -Rk.ni x^. (12) 

For the sliding interaction, the forces and toques are given by 

,(r2n2-nni)-nc 



\\Xy-X2\\ 

And for the twisting interaction it applies 

Ft = 0, 
Ml - -k,(f> . 

The constants A;,, k^, and k^ are given by 

4Fe 

R ' 

^s = 8aoG* , 
t 16 3 

^ T " ■ 
2.1.6. Inelastic interaction 



kr 



(13) 
(14) 



(15) 
(16) 



(17) 
(18) 
(19) 



As already mentioned before, inelastic motion sets in when the 
displacements exceed a critical displacement ^crit, (cm, or ^crit- 
Physically this coiTesponds to energy being dissipated upon re- 
location of the contact area. The critical sliding and twisting dis- 
place ments have been derived theoretically (Dominik & TielensL 
119961) as 



fecrit — 



2-v 

I 

'l6n' 



flO, 



The val ue of the critical roUi r ig disp lacement ^crit is still debated. 
At first, 'Domi nik & Tielerisl (Il997h assumed ^cm close to inner- 
atomic distances and chose fcrit - 0.2n m. However, subsequent 
laboratory experiments by iHeim et aP d 19991) indicate a much 
higher value of ^ciit = 3.2nm for sph erical Si02 grains of 1.9jum 
in diameter. In this work we follow iPaszun & Dominilg (l2008h 
and set ^crit = 2nm for 1 .2jjm sized grains. 

When a displacement exceeds its critical value it will be re- 
stored to the elastic limit. If for instance the rolling displace- 
ment exceeds its critical value ||^|| > ^crit, the contact pointers 
will be corrected ni — > m^ n2 — » Mj and thus ^ ^> ^ such that 
ll^'^ll = ^crit- This modification of the contact pointers reflects a 
change of the corresponding potential energy. Therefore, we can 
keep track of the amount of dissipated energy. For fu rther details 
on how the inelastic motion is applied we refer to IWada et al] 
(.2007.) . 



2.2. Particle-Wall interaction 

To study the compression of a sample aggregate it is necessary 
to constrain the motion of the monomers using suitable bound- 
ary conditions. In the corresponding experiments the sample 
has been put in a solid box with a movable piston on top (see 
Guttler et al., 2009, Fig. 2). In our simulations we model the ex- 
perimental setup by putting the sample aggregate into a box of 
fixed walls. During the simulation the top wall is moving down- 
wards with constant speed in order to compress the sample (see 
Fig.EJ. 

We assume that the monomers may interact with a wall 
in a similar way as they interact with other grains. In accor- 
dance with the particle-particle interaction we derive the cor- 
res ponding forces and torques following the approach presented 
bv' Wada et al.l (l2007l) . For this purpose we assume that the wall 
can be described as a very huge particle in the limit rwaii — > °°- 



2.2.1. Compression 

The force of the particle-wall interaction in normal direction is 
very similar to the case of the particle-particle interaction. The 
compression length 5 is given by 

6 - r - d , 

where r is the radius of the monomer and d denotes the distance 
between the surface of the wall and the center of the grain. Given 
a point p located on the surface of the wall and the surface nor- 
mal iIk, we can easily obtain d by 

<i= \(x-p)-n„\ , 

where x denotes the position of the particle. The force in normal 
direction can then be calculated using Eq. ([3]). Note that here the 
reduced radius R is different to the case of particle-particle inter- 
action. In the limit of r„aii -^ co the reduced radius R equals the 
radius r of a monomer: 

„ ,. r r„aii 

R = lim r . 

'■waii-»'» r + r„aii 



2.2.2. Rolling 

Keeping in mind that R - r, the torque Mr acting on the particle 
caused by rolling along the surface of the wall is given by 



Mr = ^r.wall?- «1 X «„ 



(20) 
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Fig. 2. Setup of the numerical simulations to obtain material parameter for porous dust agglomerates. Left: Omnidirectional com- 
pression, while the top wall is moving downwards at constant speed, the sample is enclosed in a box with fixed walls on all sides. 
Right: Unidirectional compression, as the sample is getting compressed between two walls, particles can leave the initial volume 
to the sides. Particles colored red are actually in contact with the walls. 



where ^r.waii is equivalent to the rolling constant k^ given in 
Eq. (fTTb taking the different reduced radius of the particle-wall 
interaction into account. «„ denotes the surface normal of the 
wall. Note that it is important on which side of the wall the par- 
ticle is located with respect to the direction of «„. In the sim- 
ulation we must either ensure that the particles remain on the 
positively oriented side all the time or check on which side of 
the wall the particle is located and correct the sign of Eq. ( |20| ) if 
necessary. 

To model the motion of a particle which is rolling inelasti- 
cally over a wall we use a similar approach as for the inelastic 
particle-particle interaction. Taking r2 — > oo into account, the 
correction of the contact pointers is then given by 



111 = «i - -A^ , 



r 



(21) 

Mj = ni-n^, (22) 

where a is a correction factor derived in detail in IWada et aP 
( l2007h and A^ is given by 

Af = f(l-|;)^ (23) 

The contact pointer «2 of the "wall particle" is equivalent to the 
normal vector n^ of the wall and is not modified during the in- 
elastic rolling motion. 



2.2.3. Sliding 

To describe the sliding motion of a particle on a wall it is not suit- 
able to start with Eq. (|9]l and assume that Mc — > «,, and M2 — * n„. 
Therefore we choose a different approach that takes into account 
how far the contact has slided from the location where it has 
initially formed. 

Let p denote the position, where the surface of the particle 
touched the wall when estabhshing the contact. For any later 
time, the center of the contact area is given hy x + riii, where x 
is the current position of the particle. We define 

(Q^x + rtix-p . 



The sliding displacement is then given by 

^ = ^0 - (fo • «w) "w ■ 



(24) 



To model the inelastic wall sliding we modify the initial cen- 
ter of the contact area p. If ||^|| > ^crit, we apply the correction 



""'"^'■-pK- 



2.2.4. Twisting 



(25) 



The torque caused by the twisting motion is calculated in the 
same way as if two particles are in contact. Starting with the 
twisting displacement given in equation ( fTOb we obtain 



J to 



(26) 



under the assumption that the wall does not rotate. The torque is 
then given by 



Mt.wall - -^t,wall0wall 



(27) 



where fctwaii = K- Here we assume that particles may twist elas- 
tically around the same angle for both the particle-particle and 
particle-wall interaction. According to test simulations where we 
measured the relative importance of the different wall interaction 
dissipation channels, inelastic wall twisting is of only minor im- 
portance. 



3. Setup 

In order to calibrate the model and to compare simulations in de- 
tail with experimental results we focus here on two different nu- 
merical experiments that follow closely the experimental setup. 
Specifically, we will deal with agglomerates enclosed in a box 
and aggregates confined between two plates as depicted in Fig.|2] 
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3.1. Sample generation 

In accordance with the i nentioned laboratory experiments by 
iBlum & Schrapletj (|2004 : iGiittler et al] (|2009^, our samples are 
produced by random ballistic deposition (RBD), which means 
that single grains are successively poured down on the existing 
aggregate impacting from the same direction. In order to prevent 
any restructuring upon impact the impact velocity of a monomer 
hitting the sample is kept very low. The resulting samples feature 
a filling factor between 0.12-0.15. Filling factors of 0.12-0.14 
result from the fluffier surface and are therefore only observed 
for small samples. As the size of the samples increases this sur- 
face effect becomes negligible and the filling factor converges 
to 0.15 which agrees well with the dust cakes used in the corre- 
sponding laboratory experiments and with numerical studies by 
IWatsonetall(ll99 



3 



3.2. Measurements 

3.2.1. Pressure 

In the numerical (and laboratory) experiments a box-shaped 
sample is enclosed between walls that constrain the motion of 
the grains (see Fig.|2]l. Then, the top wall is being moved down- 
wards at a constant speed until the filling factor exceeds a cer- 
tain threshold 0crit- Typically, 0crit is set to 0.7. As the top wall is 
moving downwards the volume of the sample decreases and it is 
increasingly more compressed. The total force F^ acting upon 
the top wall is calculated by summarizing the forces F, exerted 
on the wall by grains that are currently in contact with it, where 
only the component in normal direction n^ to the wall is taken 
into account 



/^ w — "w 






The pressure P is then given by 

A 

where A denotes the base size of the box. 

If there are only a few particles in contact with the wall, 
Fh, may change considerably from one integration step to the 
next due to the normal oscillations of the particle-wall contact 
(see Sect. l2.1.4l i. Since these vibrations occur on a timescale 
of nanoseconds whereas the compression timescale is typically 
orders of magnitudes higher, it is reasonable to average over 
several integration steps (covering a few normal particle os- 
cillations) to reduce the noise in the pressure determination. 
Typically, we averaged over 100 integration steps in this work. 



3.2.2. Filling factor 

The volume filling factor (p is defined as 



<f>. 



V, 



mat 

n7' 



(28) 



where Vmat = AjTinr'N denotes the volume of A^ particles of 
radius r and Vtot is the volume that the aggregate currently occu- 
pies. Calculating the filling factor is trivial in the case of omni- 
directional compression (see Fig.|2] left panel) as Vtot = Ah for a 
box with base size A and current height h. 

However, in the uni-directional case (Fig.|2] right panel) 
there are no side walls containing the aggregate, and particles 



will leave the initial volume. They flow to the sides as the top 
wall is moving downwards. This complicates the determination 
of Vtot in the numerical as well as experimental setup. In the 
following we assume that the volume the aggregate is currently 
occupying is given by its projected cross section Aproj and the 
cuiTent height h of the aggregate, which equals the distance be- 
tween the top and the bottom wall. Ap,oj is obtained by projecting 
the aggregate in the plane of the top/bottom wall. 



3.3. Previous work 

A first attempt to determine the compressive str ength numer- 
ically has been presented by iPaszun & Dominilg (12008.) . With 
the exception of the damping of the normal interaction (see 
Sects. |2Tj1 and 12.1.4b they used the same particle interaction 
model that we use here. However, instead of flat walls they used 
two huge particles with radii much bigger than the dust grains to 
model the boundary conditions. The sample was put between the 
"wall particles" and the upper particle was being moved down- 
wards at constant speed while the force acting on this "wall par- 
ticle" was stored for later analysis. Since the motion of the grains 
has only been constrained in the vertical direction, grains could 
dodge to the sides as the pressure increased. Thus, an increase 
of the initial cross section of the aggregate was observed during 
the compression that led to a significant uncertainty in the deter- 
mination of the volume the aggregate occupied at a certain point 
of time. Since the number of monomers was also limited to very 
smaU numbers (a; 300) it is questionable if the results hold in the 

co ntinuum limit. 

iPaszun & Dominilg (2008) confined themselves to the case 
of unidirectional compression. To our knowledge, the case of 
omnidirectional compression has not been simulated so far for 
this material. The compressive strength was only determined for 
a compression speed of 0.05 m/s. A possible dependency of the 
compression behaviour on the speed of the compression has not 
been studied before. 



4. Numerical method 

To integrate the equations of motion we use a second order, sym- 
plectic Leap-Frog scheme. The main reason behind this choice 
lies in the fact that the forces and torques have to be deter- 
mined only once during each integration step. Likewise, the evo- 
lution of the contact pointers and the twisting displacement (see 
Eq. ( [Tol l) is calculated with second order accuracy. 

The timestep is limited by the oscillations in direction of the 
contact normal caused by the normal force (see Sect. l2.1.4l i as 
well as the oscillations in the tangential plane of the contact 
caused by the sliding force. As already mentioned the normal 
oscillation period is of the order of 10 ns. 

The period of the oscillati ons in the ta ngential plane can be 
obtained in the following wav: IWada et al.k2007 ) derive the slid- 
ing force from the corresponding sliding potential L'^s 

Us = ^kMf ■ (29) 

We can get an estimate of the oscillation period T by 

T^ — ^lnJ-, (30) 

(Jj \ fcs 

where co denotes the corresponding angular frequency of the os- 
cillation and m is the mass of a monomer. For the material pa- 
rameters given in Tab. [1] we obtain T - 12.8 ns. Applying the 
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1 0^ 1 0^ 

pressure [Po] 



Fig. 3. The compressive strength (filHng factor 4> versus pres- 
sure p) as obtained from experiments. The dark shaded region 
with the soUd line fit ref er to the omn i dire ctional compression 
experiment performed bv lGiittler et alJ (l2009i) whereas the light 
shaded region and the da shed line fit were obtaine d for unidi- 
rectional compression by Blum & Schrapler (2004). The small 
image on the top left depicts the experimental setup of the 
omni-directional co mpression experiment. (Figure taken from 
lGuttleretal.[l2009h 



sliding modifier nis - 2.5 (see Sect. l5.lb the tangential oscilla- 
tion period decreases to 7" = 7.83 ns, which limits our timestep 
to 0.3 ns. 



5. Results 

Calibration experiments using a continuum SPH-code indicate 
that the compressive strength of a porous du st aggregate depends 
on how fast the compac tion takes place JGiittler et al.L 120091; 
iGere tshauser et all 1201 0^. So far, in laboratory experiments the 
compressive strength could only be determined for a slow quasi- 
static compression process, where the compressed aggregate had 
been given sufficient time for relax ation (in the following re- 
ferred to as static compression) (Gutt ler et al.L 12009). 

The static compression provides us with the possibility to 
check how well our model is able to describe the compression 
behaviour of porous dust aggregates. In the first step we will 
therefore use the case of the omnidirectional static compression 
to calibrate our molecular dynamics model. Afterwards we will 
increase the speed of the top wall. At a sufficiently high speed 
the relaxation of the aggregate will not be possible any longer 
(from now on referred to as dynamic compression). 



5. 1 . Calibration of our model 

To model the quasi-static compression the speed v^aii at which 
the top wall is moving downwards should be as low as possi- 
ble. As the time required to reach a certain filling factor is in- 
versely proportional to Vwaii, the runtime of the simulat ion con- 
stitutes a lower limit of Vw^n. lPaszun & DominikI (l2008h consid- 
ered v„aii = 5 cm/s to be low enough to model static compres- 
sion. However, we observe that the curves still change when us- 
ing even lower velocities (see Fig.|6]below). To model the case of 
static compression we set y^aii - 1 cm/s. To ensure it is a reason- 
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Fig. 4. Compressive strength for omnidirectional compression of 
cubical aggregates. The result of the unaltered interaction model 
is compared to our improved model. The black hne repres ents a 
fit to experimental data obtained from lGiittler et al.l (l2009h . 



able choice we checked lower velocities down to v„aii - 0.2 m/s 
but observed only a tiny deviation of the resulting curves. 

Depending on the number of particles we use, the diameter 
of our sample aggregates is roughly 40 - 60 //m which is about 
xlO"* times smaller than the samples used in the laboratory ex- 
periments. While the sample is getting compacted the particles 
on the edges of the sample must overcome the sliding resistance 
of the side wall. Owing to the small diameter of the sample this 
has a significant influence on the resulting force on the top wall. 
To mitigate this effect we reduce the strength of the rolling, slid- 
ing, and twisting interaction between particles and the side walls 
by a factor of 1000. 

The results of the corresponding laboratory experiment are 
shown in Fig. [3] The solid black curve is a good fit to the omni- 
directional experimental data (see Fig.|4]and dark shaded area in 
Fig.O and is given by 



cf>(P) = 02 - 



»2 -<Pl 



exp( 



loglpf-logloAT 



)-' 



(31) 



Using 01 = 0.15 and <p2 = 0.58 we obtain p^ = 16.667 kPa and 
A - 0.562. The black curve shown in Fig.|4]depicts this fit. 

The results of the unidire ctional experiments , also d isplayed 
in Fig. [3] have been fitted bv iBlum & Schraplej (1200 4) using a 
similar curve with <pi = 0.15, <p2 = 0.33, pm = 5.6 kPa, and 
A = 0.33. 



5.2. Omnidirectional compression 

5.2.1 . Quasi-static case 

In Fig.m we compare the experimental fitting curve to results 
from our simulations of box-shaped samples featuring approxi- 
mately 11000 particles and an edge length of about 40yum. All 
curves have been obtained from averaging the results of six inde- 
pendent runs with samples having statistically equal bulk prop- 
erties. 

Obviously, simulations using the Dominik and Tielens 
model do not reproduce the experimental data very well. For low 
pressures the blue dashed curve in Fig.|4]lies significantly above 
the solid black one, i.e. applying the same pressure the sample 
has been compressed to a higher filling factor in the simulations 
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than compared to the experiments. In order to make it harder to 
compress an aggregate we tried to increase the stiffness of the 
aggregates by modifying the particle interaction. 

To our knowledge the equations describing the rolling and 
sliding interaction have not been experimentally tested yet, 
whereas the pull-off force of the normal i nteraction has bee n 
measured using atomic force microscopy dHeim et al.L 1 19991) . 
Thus, we vary the strength of the rolling and sliding interaction. 
For this purpose we simply multiply the constants k^ and k^ (see 
Eqns.fTTjandfTSTl with correction factors that we further refer to 
as rolling/sliding modifiers m^ and mj. This constitutes a straight- 
forward approach to increase the stiffness of monomer chains. 
In fact, we also modified the strength of the twisting interaction 
but found that it had very little to no impact on the compressive 
curve. Therefore we do not alter twisting in this work. 

Choosing m, - 8 and m^ = 2.5, we obtain the red-dotted 
curve in Fig.|4l All in all, our modified interaction model is able 
to reproduce experimental results much better than the original 
version. In particular, for pressures below lOOkPa we observe a 
very good agreement with experimental data. However, we ob- 
serve a deviation for pressures above 300 kPa. In our simulations 
a pressure of more than 1 MPa is required to further compress 
aggregates when a filling factor of about (D = 0.52 is reached. 

In the quasi-static case the aggregate is given sufficient time 
to restructure and counteract the pressure exerted on the walls. 
Thus, we expect the filling factor increases homogeneously in 
the sample. In Fig. |5] the vertical profile of the filling factor is 
plotted for different stages of the compression process. To de- 
termine the filling factor profile, the sample is split vertically 
into equidistant intervals with the length of one particle diame- 
ter. Then, the average filling factor is calculated for each inter- 
val. Note that the fractal structure resulting in a filling factor of 
(f) - 0.15 in the bulk part of RBD-generated aggregates is not 
present at the bottom of the aggregate. Therefore the filling fac- 
tor there exceeds the average value of </> = 0.15. 

During the compression process several snapshots of the fill- 
ing factor profile have been taken. As expected the filling factor 
increases almost homogeneously for slow compression speeds. 
Keep in mind that it requires higher pressure to compact particles 
close to a wall since neighbouring particles have to be pushed 
away. Therefore the filUng factor of the particle layers close to 
the top or bottom wall is lower compared to the rest of the ag- 
gregate for highly compacted aggregates. This effect causes the 
crescent shaped curves observed for highly compacted aggre- 
gates in Fig.|5] 

The coiTesponding compressive strength curve for these low 
velocities is shown in Fig.|6] For compression speeds of 5 cm/s 
and below the differences to the quasi-static case remain small. 
For higher velocities we observe an incre asing deviation from 
the quasi-static curve. Thus, the results from lGiittler et al.l (l2009l) 
cannot be applied for higher velocities. 

5.2.2. Dynamic case 

In a second step, we studied the influence of large compres- 
sion speeds. If the compression speed exceeds around 1 m/s the 
compression behaviour changes considerably. The required sim- 
ulation time is inversely proportional to the compression speed. 
Therefore we use a higher number of particles for compres- 
sion experiments with wall speeds above 1 m/s. The box-shaped 
samples are composed of about 40000 particles and feature a 
base length and height of x 60jum. Compared to the quasi-static 
case the shape of the curves changes drastically, see Figs. Q and 
[8] Instead of a smooth transition, three distinguished regimes 
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Fig. 5. Snapshots of the vertical profile of the aggregate's filling 
factor averaged horizontally. Shown are results at 5 evolutionary 
times using 3 different speeds. The left (solid) curves represent 
the initial state at f = 0. As the top wall moves slowly down- 
wards the aggregate is compacted almost homogeneously as is 
indicated by the nearly vertical curves. At the top and bottom of 
the box wall effects produce a slight deviation. 
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Fig. 6. Compressive strength for low velocity omnidirectional 
compression. 



emerge: The filling factor does not increase until a certain criti- 
cal pressure is reached. Then, only a small additional pressure is 
required to compact the aggregate. When the aggregate is close 
to its final compaction the pressure again increases sharply. 

This can be easily explained by looking at Fig.|9] When the 
compression speed exceeds a value of 1 m/s the aggregate is 
compacted inhomogeneously. The compression occurs too fast 
to allow the propagation of the top pressure through the entire 
sample. We clearly see the emergence of a very dense layer right 
beneath the top wall. This compact layer propagates downwards 
at the speed of the wall similar to a snowplough clearing freshly 
fallen snow. While pushing the dense layer downwards the pres- 
sure remains constant. After it reaches the bottom of the sample 
the pressure required to compress the sample a little bit further 
increases drastically. The sharp spikes shown in Figs.|7]and[8]re- 
sult from the highly compacted layer reaching the bottom of the 
sample. The density wave is reflected from the bottom causing 
heavy fluctuations of the pressure on the top and bottom wall. 
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Fig. 7. Dynamic omnidirectional compression with v„aii between 
1 - 5m/s. 




Pressure [Pa] 



Fig. 8. Dynamic omnidirectional compression with v„aii between 
7 - 15m/s. 



By comparing the filling factor profile during the compres- 
sion to the initial one we can determine at which speed Vp the 
compaction is propagating downwards trough the sample. For 
this purpose we measure the height where the initial and cur- 
rent filling factor profile deviate from each other and use the 
time that has passed since the start of the simulation to calcu- 
late the speed. Averaging over six different samples we obtain 
vp = 6.98 ±0.16 m/s. 

To provide continuum-simulations with a simple recipe for 
the dynamic compressive strength we performed simulations us- 
ing compression speeds up to 25 m/s. A few examples are shown 
in Figs.|2]and[8] For every compression speed we determine a fit 
curve similar to Eq. (l3Tl i. where p,„ and A serve as fitting param- 
eters. Thus, we obtain values of pm and A for different compres- 
sion speeds. In the last step we determine for each parameter an 
analytic approximation describing the dependency on the com- 
pression speed. 

We observe different behaviour between the the quasi-static 
case for low velocities and the dynamic case for higher veloci- 
ties. In the beginning, p,„ decreases which means that the sam- 
ple can be compressed more easily. This can be explained by the 
fact, that the aggregate is given less time to restructure and coun- 
teract the external pressure exerted on it by the wall. However, 
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Fig. 9. Snapshots of the vertical profile of the filling factor. As 
the top wall moves downwards rapidly the compaction of the 
lower parts of the aggregate is lagging behind. The color indi- 
cates the compression speed whereas the line type indicates the 
position of the top wall. The solid lines show the filling factor 
profile of the initial uncompressed sample. 



this effect will be reversed when the compression speed exceeds 
a critical value of Vcrit ~ 0.9 m/s. In the dynamic regime, it gets 
considerably harder to compress the sample with increasing ve- 
locity of the wall. Therefore it is helpful to distinguish between 
the two regimes. 

In Fig.[10] the dependence of /?„, on the compression speed 
v is shown for values of v < 1 m/s. Using the ansatz Pm{v) - 
av^ + bv + c we obtain 



p^{v) = (18.296 v^ - 33.663 v + 16.667) kPa , 



(32) 



where the compression speed v is given in m/s. 

In Fig. [TT] the dependence of p,n on the compression speed 
V is shown for values of v > 1 m/s. To find a simple analytic 
approximation we choose a power law of the form Pmi}') - av'' + 
c and we obtain 



1.93 



p„iv) = (1.340 v'-^" -H 0.307) kPa 



(33) 



The fitted parabola describes the data points well. As the expo- 
nent of 1 .93 is close to 2, we determine a second fit where the 
exponent was set to 2 and get 

pjv) = ( 1 .087 v^ + 0.560) kPa . (34) 

Similarly, for the parameter A (see Fig.[T2li we obtain 

A(v) = (v + 1.598)"'-''^^ + 0.170 . (35) 



5.3. Unidirectional compression 

Additionally, we simulated the unidirectional compression of 
cylindrical samples of different sizes using the non modified 
model. The results are shown in Fig.[T3]where again each curve 
is obtained by averaging the results f or six different samples of 
equal size. To compare our results to lPaszun & DominikI (l2008h 
the speed of the top wall was set to Vwaii = 5 cm/s. Apparently 
there is a noticeable discrepancy between our simulations and 
the laboratory results. As in the case of omnidirectional com- 
pression, the pressure required to reach a certain filhng factor is 
significantly lower in our simulations. 
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Fig. 11. Dependence of the fit parameter p„ of Eq. (ISTT i on the 
compression speed v in the dynamic regime of v > 1 m/s. 



As we can see in Fig. [13] the deviation for pressures above 
lOkPa becomes more apparent if we increase the size of the 
samples. As Paszun & Dominik (2008) compressed very small 
samples using only about 300 particles this may be the reason 
why their results showed better agreement with laboratory re- 
sults for higher pressures. However, their compressive strength 
curve was also shifted in the same direction as in this work. 

Afterwards we tested the modified model with the same 
m, - 8 and m^ - 2.5 as found above. As it is shown in Fig. [14] 
(red dotted curve), the modified model agrees very well with 
the experimental results for pressures below lO'^Pa. However, 
we still end up with considerably higher filling factors. A pos- 
sible explanation is given by the fact that the size of the sam- 
ples used in the laboratory experiments was of the order of cen- 
timeters and thus about 2000 times larger than our samples. 
[Blum & Schriipler (2004) measured that the projected cross sec- 
tion increased by a factor of 1 .6 during the compression process. 
To reach similar pressures in our simulations the sample has to 
be compressed until its height is only about two times the di- 
ameter of a single monomer where the cross section increased 
roughly by a factor of 4. Due to the larger diameter of the labo- 
ratory samples it is harder for monomers in the center to flow in 
the outward direction. 
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Fig. 13. Compressive strength for unidirectional compression 
of cylindrical aggregates of different size using the original 
Dominik & Thielens model. The black line has been obtained 
from experiments by Blum & Schrapler (2004) . 
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6. Conclusions 

We have performed molecular dynamics simulations to study 
the compressive strength of dust agglomerates which plays an 
important role in determining the outcome of mutual collisions. 
Using a special setup for the simulations we were able to com- 
pare our results in detail to the outcome of of dedicated labora- 
tory experiments. 

Our simulations using the frequently applied interaction 
model by Dominik & Tielens (1997) indicate that real aggre- 
gates composed of micron sized silicate grains feature a greater 
stiffness. Since the primary bulk properties of material used for 
the individual monomers are known experimentally very well 
(see Table 1), we decided to vary the force constants (A:, and kg) 
for rolling and sliding. Indeed, the higher stiffness can be accom- 
modated by an increase of m,- - 8 and m, = 2.5 (for kr and ks, 
respectively) in comparison to the quoted values in Eqns. (fTTT i 
and ( fTST l. After modifying the interaction model as described 
in Sect. 15. II we have been able to reproduce the experimental 
results much better, and found very good agreement for both, 
unidirectional and omnidirectional compression. This work re- 
veals the importance of the rolling and sliding interaction for the 
restructuring of aggregates. As these interactions currently lack 
experimental testing we feel it desirable to study in particular the 
rolhng and sliding of micron sized grains in laboratory experi- 
ments. 

We have also studied the influence of the wall speed on the 
compression behaviour. If an aggregate is compressed slowly the 
filling factor increases homogeneously and the pressure needed 
to further compact the aggregate increases with increasing fill- 
ing factor. For higher compression velocities a compacted layer 
emerges underneath the moving wall, similar to the shovel of a 
snow plow when pushing away snow. Once this layer has formed 
the pressure remains nearly constant until the layer has reached 
the bottom of the sample. The transition from the static towards 
the dynamic case occurs at compression speeds of the order of 
1 m/s. Since impact velocities of typical collisions of planetesi- 
mals lie within the range of 1 - 10 m/s the dynamic compression 
behaviour must be taken into account when simulating such col- 
lisions. 

To determine the impact of the rescaling of the rolling and 
sliding forces on the very early phases in the planetesimal forma- 
tion process, we plan to perform detailed collision simulations 
for a wide set of collision parameter. This will allow us to calcu- 
late ab-initio the division between bouncing, sticking and frag- 
mentation. This has recently been under experimental and th eo- 
retical scrutiny dZsom et aI.L l20ld : iGeretshauser et all 1201 Ih . It 
is hard to estimate the consequences of the new (stiffer) parame- 
ter set on the outcome of agglomerate collisions. We suspect that 
this can possibly lead to an enhanced sticking as more energy 
can be stored in the system before it breaks. Using a much larger 
particle number than previously will allow us to determine more 
accurately important bulk parameter such as the sound speed. 
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